#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>

//	for (i = 0; i < n*dim; i += 1) printf("%0.f ", X[i]); printf("\n\n");
//	t1 = omp_get_wtime();
//	t2 = omp_get_wtime();
//	printf("Step 1 - %g secs\n", t2-t1);

// blockBits is the number of bits for a block, presumable fits in cache.
void hybridScan(double *X, int n, int blockBits, int dim) {
	int i, j, k, d, logn, step, limit, mask, bit, st, maxIdx, srcIdx, targetIdx, bitdim;
	double t1, t2;
	double *srcPtr, *targetPtr;
	logn = (int)ceil(log2(n));
	
	step = 1 << blockBits;  // shift the 1 into the blockBits+1 position
	//for(j = 0; j < n * dim; j++) printf("%.0f ", X[j]); printf("\n");
	
	//t1 = omp_get_wtime(); 

	  // Split n into blocks of 2^blockBits
	  // Do sequential scan on each block
#pragma omp parallel for shared(X) firstprivate(dim, step, n) private(i, st, j, maxIdx, srcPtr, targetPtr) schedule(guided)
	for (i = 0; i < n ; i += step ) {
		st = (i + step) > n  ? (n - i) : step ;
		maxIdx = (i + st) * dim;
		for (j = (i + 1) * dim; j < maxIdx; j +=dim) {
			targetPtr = X + j;  // no need to multiply by sizeof(double) since the pointer is double*
			srcPtr = X + j - dim;
			for (k = 0; k < dim; k++)
				targetPtr[k] += srcPtr[k];
		}
	}
//	t2 = omp_get_wtime(); 
//	printf(" forward time = %g\n", t2 - t1);
	//for(j = 0; j < n*dim; j++) printf("%.0f ", X[j]); printf("\n");

//	t1 = omp_get_wtime(); 
	  // next do a parallel scan on only the "splitter" elements
	  // step = 1 << blockBits;
	mask = step - 1;
	bit = step;

	for (d = blockBits + 1; d <= logn; d++) {
		mask |= bit;  // increase mask
		step = bit << 1;
		  // has to have at least mask to enter if.
#pragma omp parallel for shared(X) firstprivate(mask, n, step, dim, bit) \
	private(i, j, targetIdx, srcIdx) schedule(guided) 
		for (i = mask; i < n; i += step) {
			targetIdx = i*dim;
			srcIdx = (i^bit) * dim;
			for (j = 0; j < dim; j++)	
				X[targetIdx + j] += X[srcIdx + j];	// now do the math
		}
		bit = step;
	}

	mask = (1 << logn) - 1;
	bit = 1 << (logn-1);

	// skip top bit.  each run updates 1 below.	
	// since each run udpates 1 lower in bit level,  
	// d == blockBits + 1 would update d = blockBits, 
	// which would be used to udpate the other entries in blockBits.
	for (d = logn-1; d > blockBits; d--) {  
		mask ^= bit;  // keep lower bits only
		step = bit;
		bit >>= 1;   // use to find the half step merge target.
		limit = n - bit;
		
		// update i, i is the paired target.  j has the computed sum that
		// needs to be added to i. e.g. i = 3, and j = 2.  this updates value
		// of i.
		
		// has to have at least mask, and can't go over n-bit.
		bitdim = bit * dim;
#pragma omp parallel for shared(X) firstprivate(mask, limit, step, dim, bitdim) \
		private(i, j, srcIdx, targetIdx) schedule(guided) 
		for (i = mask; i < limit; i += step) {
			srcIdx = i * dim;
			targetIdx = srcIdx + bitdim;
			for (j = 0; j < dim; j++)	
				X[targetIdx + j] += X[srcIdx + j];
		}
	}

//	t2 = omp_get_wtime(); 
//	printf(" middle time = %g\n", t2 - t1);
	//for(j = 0; j < n* dim; j++) printf("%.0f ", X[j]); printf("\n");

//	t1 = omp_get_wtime(); 
	
	  // finally do blockwise sequential update at the very end.
	step = 1 << blockBits;  // reset step size to block size.	
#pragma omp parallel for shared(X) firstprivate(dim, step, n) private(i, st, j, maxIdx, srcIdx, srcPtr, targetPtr) schedule(guided)
	for (i = step ; i < n ; i += step ) {
		  // not >=, which causes the n-i branch to be taken.
		srcPtr = (double*)malloc(dim * sizeof(double));
		st = (i + step ) > n  ? (n - i) : (step - 1);   
		srcIdx = (i - 1) * dim;
		for (k = 0; k < dim; k++) {
			srcPtr[k] = X[srcIdx + k];
		}
		maxIdx = (i - 1 + st) * dim;
		for (j = i * dim; j <= maxIdx; j += dim) {
			targetPtr = X + j;
			for (k = 0; k < dim; k++)
				targetPtr[k] += srcPtr[k];
		}
		free(srcPtr);
	}
//	t2 = omp_get_wtime(); 
//	printf(" backward time = %g\n", t2 - t1);
	//for(j = 0; j < n * dim; j++) printf("%.0f ", X[j]); printf("\n");
}
